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ABSTRACT 

Windowing of design space is considered in order to 
reduce the bias errors due to low-order polynomial 
response surfaces (RS). Standard design space windowing 
(DSW) uses a region of interest by setting a requirement on 
response level and checks it by a global RS predictions 
over the design space. This approach, however, is 
vulnerable since RS modeling errors may lead to the wrong 
region to zoom on. The approach is modified by 
introducing an eigenvalue error measure based on point-to- 
point mean squared error criterion. Two examples are 
presented to demonstrate the benefit of the error-based 
DSW. 

1. INTRODUCTION 

The popularity of response surface (RS) techniques in 
design optimization studies has brought attention to ways 
of increasing the accuracy of RS approximations. 
Adequacy and accuracy of RS are mainly affected by the 
following three factors; 

• Use of finite number of data points due to cost of 
data generation 

• Noise in the data 

• Inadequacy of the fitting model 

We are mainly focused in this paper on model 
inadequacy or bias error due to use of low-order 
polynomials such as quadratic RS. 

An obvious way for reducing bias error is to use higher 
order polynomials or more complex functions in RS. 
Cubic or even higher order polynomials were applied, for 
instance, by Venter et al. [1] and by Papila (N) et al. [2]. 
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Papila and Haftka [3-4] also achieved substantial 
improvement in accuracy by using cubic polynomials in 
RS approximation for HSCT wing bending material 
weight. 

A high-order fitting model, however, requires a large 
number of data points, which is usually prohibitive in high 
dimensional problems. 

Design of experiments (DOE) offers minimum-bias 
criterion for reducing modeling (bias) error when low- 
order models are used. Venter and Haftka [5] developed 
an algorithm implementing minimum-bias based criterion, 
necessary for an irregularly shaped design space where no 
closed form solution exist for minimum-bias design. 

Decreasing bias error is also possible by reducing the 
size of the fitting region by the use of reasonable-design- 
space (RDS) approach. The approach starts by identifying 
constraints specific to the problem of interest provided that 
they are easy and inexpensive to evaluate. These 
constraints are then used to eliminate unreasonable designs 
from the original design space. For instance, simple 
geometric constraints where applicable may prevent 
combinations of design variables resulting in unreasonable 
geometry configurations [4, 6-8]. Finally, tools of DOE 
such as D-optimality select data points within the region of 
interest where response is evaluated. 

It is also possible to identify region or regions of 
interest by windowing the design space simply based on 
the observed or predicted response levels over the design 
space [2]. The windowing approach shares the sprit of 
RDS approach in that it reduces bias errors by reducing the 
size of the design domain. However, windowing needs 
data generation and a global RS beforehand unlike 
traditional RDS approach. 

In this paper, we aim to focus on improving the RS 
accuracy particularly in the regions critical to design 
optimization. We limit ourselves to quadratic RS 
approximation and concentrate on effective use of design 
space windowing (DSW) approach relying on a global RS. 
Standard DSW uses a requirement set on the response 
level and checks it by using a global RS predictions over 
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the whole design space. This approach, however, is 
vulnerable since poor accuracy in the global RS may lead 
us to zoom on the wrong region. 

At this point, we modify the approach and call it as 
error-based DSW by introducing a bias error measure. 
The error measure is based on an eigenvalue problem 
obtained by point-wise mean squared error criterion. 
The eigenvalue problem was derived by Papila and Haftka 
[9] and used as a tool to map qualitatively the RS bias error 
by the associated eigenvalues. In the proposed approach, 
regions with high eigenvalues corresponding to potential 
large bias error are excluded from consideration. 

The following section presents more detailed 
description of the mean squared error criterion and use of 
eigenvalues in the DSW approach. Section 3 describes 
example problems used to demonstrate the approach. 
Section 4 presents the results and discussion followed by 
concluding remarks in Section 5. The derivation of 
eigenvalue problem along with the background for RS 
methodology can be found in the Appendix. 

2. APPROACH 

2.1 Mean Squared Error Criterion-Eigenvalues 

An approach for estimating RS approximation bias 
errors due to fitting model inadequacy is presented in Ref. 
[9], The mean squared error predictor ( MSEP ) for an 
inadequate model is studied point-to-point yielding an 
eigenvalue problem where the maximum eigenvalue at 
each point provides a relative estimate of maximum bias 
error. With the calculation of the maximum eigenvalues 
over the design space, regions of possible high bias error 
are identified. The derivation presented in Ref. [9] can 
also be found in this paper as an Appendix. As can be seen 
from the derivation the eigenvalue error measure strongly 
depends on the DOE used, but not on the response data. 

Papila and Haftka [9] used face-centered central 
composite design (FCCD) and demonstrated the use of 
eigenvalue estimate of bias error for problems where the 
true function is a cubic while the fitting model is quadratic. 
In particular, positive correlation between the square-root 

of maximum eigenvalues and the absolute residuals 

was found for the examples of 2D polynomials that were 
studied. 

We also use FCCD as our original DOE, quadratic RS 
as our fitting model and calculate the eigenvalues as if the 
true function is a cubic. Figure 1 presents the FCCD points 

and relevant field. 

2.2 Design Space Windowing Approach 

We adopted two types of DSW approaches based on 
different reasonability conditions while windowing for the 
design region or regions of interest. For simplicity, we 
consider problems where we are interested in the high 
response regions. 


The standard DSW uses Eq. (1) as the reasonability 
condition. 

y RS\ - ymteresl 0) 

where is the prediction by the global RS (RSI) and 
y j n t erest * s response bound to define the region of 
interest. 

Error-based DSW modifies the standard DSW by 
adding another condition based on the eigenvalues [9] that 
characterize the modeling error. In order to set a 
precaution for the possible misleading information due to 
inaccuracy of the global RS (RSI), our reasonability 
conditions in error-based DSW can be written as follows. 

yRS\ ^ -Fine erest 

and (2) 

£ mean (Vv> 

where mean(^j ) is the mean over the design space. 

Design points with eigenvalues less than their mean over 
the design space are more likely to be accurately predicted 
by the global RS (RSI). 

Figure 2 presents the generalized flowcharts for the 
standard and error-based DSW approaches. The following 
descriptions give the details of our implementation in the 
flowcharts. 

Standard DSW Approach (Figure 2a): 

Step I: We start with a standard DOE, face-centered 
central composite design (FCCD) for Data set 1 and 
construct a global RS approximation (RSI). 

Step 2: We use the global RS (RSI) to predict the 
response values on design points of a fine grid netting the 
whole design space. 

Step 3: We identify design region or regions of interest 
based on the predictions in Step 2. We build a pool of 
supposedly reasonable designs by simply disregarding the 
designs violating the condition in Eq. (1). This is done for 
each identified region in case of multiple disjoint regions 
of interest. 

Step 4: Among the designs located in the pool of Step 
3, the number of data points is reduced to a desirable 
amount by using D-optimality. We include design points 
of Data set 1 where data itself is in the region of interest 
(y ^ y interest )• This reduces the number of additional 

evaluations for creating Data set 2. We then construct new 
RS (RS2) approximation to be used over the associated 
region only. (In case of multiple regions separate RS2 s 
are constructed) 

Error-based DSW Approach (Figure 2b): 

Step 1: Same as Step 1 of standard approach 

Step 2: We use the global RS (RSI) to predict the 
response values and calculate the eigenvalues on design 
points of a fine grid netting the whole design space. 
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Step 3: We identify design region or regions of interest 
combining eigenvalue information ( ) with the 

standard DSW approach by Eq. (2). 

Step 4: Same as Step 4 of standard approach except 
denoting RS in refined regions as RS3. 

3. TEST PROBLEMS 

We investigate two different test problems to help 
assess the performance of the strategies explained above: 
(1) Two dimensional preliminary supersonic turbine 
design, and (2) A quartic polynomial in two dimensions. 
The main reason for choosing 2D problems is that we want 
to visualize easily the response surface, prediction error 
and eigenvalue field distributions. 

3.1 Preliminary Supersonic Turbine Design in Two- 
Dimension 

In our previous efforts [2, 10-14], we have studied the 
preliminary and detailed design optimization of supersonic 
turbines for reusable launch vehicle (RLV). A two-variable 
version of the two-stage turbine design problem is 
considered in this study. The design variables are the 
mean diameter, D, and RPM ( 5.081 < D(in) < 15.243 and 
21,977 < RPM < 40,814 ). They are normalized as coded 
variables X\ and x 2 in (-1,1), respectively [15]. The design 
domain of the coded variables is a square as shown in 
Figure 1. 

The numerical simulations are based on the 
aerodynamic design software called Meanline Flow Path 
Generator [10, 16] that allows rapid analyses of turbine 
flow fields. Using the overall turbine and stage input, the 
Meanline code first generates a candidate turbine flow path 
and displays a plot of the elevation view. The code then 
runs a 1-D quick aerodynamic analysis, calculating gas 
conditions, velocity triangles, and required number of 
airfoils, predicted efficiency and power output. In 
turbomachinery design problems, high efficiency and low 
weight systems are sought. For this design problem, the 
compromise between these two criteria can be quantified 
by a single response that is payload of the RLV. 
Therefore, the output or response of interest from the 
Meanline code is the change in payload compared to a 
fixed baseline design (i.e. Apay). We are mostly interested 
in positive Apay designs (i.e., Apay>0). 

3.2 Quartic Polynomial in Two-Dimension 

We wanted to mimic the efficiency data of supersonic 
turbine blade shape optimization problem [10] in terms of 
the range of response values as varying between 0 and 1. 
Since we are mainly interested in high efficiency regions in 
turbomachinery designs, we will consider y>0.7 as the 
region of interest for this problem. 

A quartic function in 2-D is devised with variables 
ranging between -1 and +1 (as coded variables). The 


quartic polynomial ranges between 0 and 1 is given in Eq. 
3 

y = 0.742+ 0.000486.x, 2 + 0.000486*, x 2 
-0.242704*2 +0.0 12646* 3 +- 0.000486*, 2 * 2 (3) 

+ 0.000486*|* 2 +0.486381*| 3 *2 

where Xi and x 2 range in (-1,+ 1) as shown in Figure 1 . 


4. RESULTS AND DISCUSSION 

We assess the accuracy of RS using mainly root-mean- 
squared (rms) error calculations. 


► rms-error Predictor: 


SO, -h ) 2 


N-n h 


• Testing rms-error : 


M-N 

Z O'. -y>) 

1( M-AT 


where N , n b , M number of data set points, number of 
coefficients in RS and number of design points in the 
design space grid, respectively. 


•Testing rms-error in y £ y mXeres t : 


rs— 


where 


K is the number of the testing data in y £ y mleres( region 

excluding the data points of the associated RS. 

Face-centered central composite design (FCCD) as 
shown in Figure la is used for fitting a quadratic global RS 
(RSI). We use 21 by 21 grid over the design space for the 
evaluation and assessment of the methods are used. 
Therefore, N= 9, n b = 6 and M- 441. K is problem 
dependent as reported on the tables. 

The eigenvalue distribution ( ) is a function of 

DOE, the fitting model and assumed true model [9]. We 
use FCCD as our DOE, quadratic RS as our fitting model 
and calculate the eigenvalues as if the true function is a 

cubic. Therefore distribution of ^X G is identical for both 
problems as shown in Figure lb. 


(D Results for Turbine Design 

The Meanline results (exact function for turbine design 
problem) are shown as a contour plot in Figure 3a. Three 
quadratic RS models were studied. The details of these RS 
models are given below. 

• RS 1: Quadratic RS based on 9 design selected by 
standard FCCD (Figure la). Its prediction error 
contours are given in Figure 3b. 

• RS 2: Quadratic RS based on 9 reasonable designs 
using standard DSW condition. Four of the points 
are original FCCD points, which satisfy Apay >0. The 
other five designs are selected from the designs of 
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positive RS 1 predictions of the Apay calculated in 2 1 
by 21 grid over the design space. Design points are 
shown in Figure 4a. 

• RS 3: Quadratic RS based on 9 reasonable designs 
using error-based DSW conditions. Four of them 
are FCCD points, which satisfy Apay>0. The other 
five designs are selected from the designs where RS 1 
predicts positive Apay and eigenvalue condition 
satisfied. Design points are shown in Figure 4b. Table 
1 summarizes the statistics of the three RS of the 
turbine problem. It appears that rms-error predictor of 
RSI is conservative as the testing rms-error is less 
than the half of the predictor. Testing rms-error in 
region of interest is 63.5, slightly higher than overall 
testing rms-error. Both approximations of windowing 
approaches, RS2 and RS3, resulted in smaller testing 
rms-error in region of interest than RSI as expected. 
There is factor of about 2.8 and 4.2 in the magnitudes 
for RS2 and RS3, respectively, compared to RSI. This 
indicates that error-based DSW approach improved on 
the standard one. The substantial increase in the 
overall testing rms-error and max. error reflects the 
fact DSW approach approximations RS2 and RS3 are 
not accurate outside of the region of interest. This can 
also be observed on the RS2 and RS3 prediction error 
contours, Figure 5a and Figure 5b, respectively. 
Comparison of error contours given in Figure 3b and 
Figure 5 demonstrate the benefit from the windowing 
approaches where we are most interested. In spite of 
the fact that maximum error in region of interest by 
RS3 is higher than RS2 (Table 1), more uniform low- 
error distribution in Figure 5b compared to Figure 5a 
and reduction in testing rms-error in the same region 
shows us that we benefit from eigenvalue information 
during windowing. 

In order to check if the eigenvalues (bias error measure) 
helped us to select design points with more accurate RSI 
predictions, we compare error at selected points in Table 2 

and Table 3. The tables report the error measure yj/ l G and 

RSI prediction errors at the RS2 and RS3 design points, 
respectively. Excluding four common design points 
coming from the original FCCD, average error (average 
error of the boldface rows) decreases from 66.4 in Table 2 
to 45.7 in Table 3. This confirms the usefulness of the 
eigenvalues for this example. 

(ii) Results for Quartic Polynomial 

Contour plot for quartic example response given in 
Figure 6a shows that we can consider two separate regions 
of interest: upper-right and lower-left quadrants of the 
square design domain. This observation and the logic 
behind the windowing procedure lead us to evaluate each 
region separately. The regions, on the other hand, are not 
completely remote and disjoint. Therefore, we also want 
to see the effect of employing the DSW approaches on the 


connected regions as a single region of interest around one 
diagonal of the domain. We first present the single 
continuous region consideration; windowing as single 
region . Then present results of windowing as multiple 
regions . 

(a) Windowing as single region 

Four quadratic RS models were studied for the overall 
design domain. The details of these RS models are given 
below. 

• RS 1: Quadratic RS based on 9 designs selected by 
standard FCCD (Figure 1). Its prediction error 
contours are given in Figure 6b 

• RS 2: Quadratic RS based on 9 reasonable designs 
using standard DSW condition. Four of the data 
points are FCCD points, which satisfy y>0J. The 
other five designs are selected from the designs of 
positive RS 1 predictions of the>> calculated in 21 by 
21 grid over the design space. Design points and error 
contours are shown in Figure 7a and Figure 8a, 
respectively. 

• RS 3: Quadratic RS based on 9 reasonable designs 
using error-based DSW conditions. Four of the data 
points are FCCD points, which satisfy y>0.7. The 
other five designs are selected from the designs where 
RS 1 predicts y higher than 0.7 and eigenvalue 
condition satisfied. Design points and error contours 
are shown in Figure 7b and Figure 8b, respesctively. 

• RS 4: Quadratic RS based on RS 1 and RS 3 design 
points (13 designs in total). Design points are shown 
in Figure 9a. 

Table 4 summarizes the statistics of the four RS. 
Unlike turbine problem, rms-error predictor of RSI (equal 
to zero) suggesting perfect fit gives a sense of security that 
is proved to be wrong by the nonzero testing rms-error 
Testing rms-error in region of interest is even higher and 
equal to 0.051. Both approximations of DSW approaches, 
RS2 and RS3, resulted in higher testing rms-error in region 
of interest than RSI. They are 0.073 and 0.063, 
respectively. This indicates that error-based DSW 

approach improved on the standard one for the single 
windowing, but neither DSW helped to increase accuracy 
compared to RSI with single windowing. Figure 8a and b 
also compares visually the two DSW approaches in terms 
of error distribution. It shows that error-based DSW 
resulted smaller errors at remote locations of the quadrants 
forming the region of interest. 

Table 5 and Table 6 show error measure and 

RSI prediction errors at the RS2 and RS3 design points, 
respectively. Excluding five common design points 

coming from the original FCCD average error (average 
error of the boldface rows) decreases from 0.0916 in Table 
2 to 0.0868 in Table 3. 

The best performance among the four RS was obtained 
by RS4 (Figure 9). Although its rms-error predictor is the 
largest we obtained smallest overall testing rms-error , 
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testing rms-error and maximum error in the region of 
interest. In other words, adding new data points selected 
by DSW to the original FCCD improved the accuracy. 

(b) Windowing as multiple regions 

We report results only on the upper-right quadrant 
since the results of lower-left quadrant are mirror image of 
the upper right. For the statistics and design points we 
consider 11 by 11 grid points in the quadrant. Two RS 
were constructed as described below 

• RS 5: Quadratic RS based on 9 reasonable designs 

using standard DSW condition. Three of them are 
FCCD points, which satisfy y>0. 7. The rest is 

selected from the designs of RS 1 predictions 
exceeding 0.7. Design points for RS 5 are shown 
Figure 10a. 

• RS 6: Quadratic RS based on 9 reasonable designs 
using error-based DSW condition. Three of them 
are FCCD points, which satisfy y>0J . The rest of is 
selected from the designs where RS 1 predicts y higher 
than 0.7 and eigenvalue condition satisfied. Design 
points for RS 5 are shown Figure 10b. 

Table 7 summarizes the statistics of the two RSs. 
Overall testing rms-error and testing rms-error in region of 
interest is lower by standard DSW than error-based DSW 
(overall: 0.042 and 0.046, respectively, and in the region of 
interest: 0.026 and 0.030, respectively). In other words, 
standard DSW did a better job compared to error-based 
approach in this example. One possible reason is that 
introducing the error-based condition may increase 
irregularity of the domain of interest. We also calculated 
the RSI testing rms-error in region of interest over this 
particular quadrant as 0.061 that shows both DSW 
approaches helped us to lower by half. Figure 1 1 presents 
the error contour for RS5 and RS6. The location of high 
errors seemed to be shifted by the error-based DSW. 

Table 8 and Table 9 present eigenvalue based error 

measure and RSI prediction errors at the RS5 and 

RS6 design points, respectively. Although we did not see 
benefit from the error-based DSW in terms of statistics, 
these tables show that eigenvalues condition in fact 
selected points with lower errors. Excluding five common 
design points average error (average error of the boldface 
rows) decreases from 0.0943 to to 0.0658. 

5. CONCLUDING REMARKS 

In this paper we focused on the model inadequacy or 
bias error due to use of low-order polynomials such as 
quadratic RS and increased accuracy by design space 
windowing (DSW). We set a requirement on the response 
based on a global RS and zoom with a local RS on that 
region. Since the modeling errors may lead us to the 
wrong region to zoom on, we integrated an eigenvalue 
error measure into the procedure and called it error-based 
DSW. Two examples were studied to demonstrate the 
benefit from error-based DSW: 


• In two-dimensional two-stage turbine problem one 
region of interest was identified. Both DSW 
approaches improved on the global RS accuracy. 
Statistics showed that accuracy obtained after error- 
based DSW increased in the region of interest 
compared to standard DSW 

• In quartic polynomial example two regions were 
identified. In each region DSW approaches offered 
substantial improvement in accuracy compared to 
global RS. The error-based DSW, however, did not 
bring improvement over the standard approach 
possibly due to increased irregularity of the region of 
interest caused by the eigenvalue condition. 

• For both examples, the average errors at the data 
points are lower after the selection by error-based 
DSW. This is an indication that modified approach 
employing eigenvalues disregarded design points 
where reasonability based on RSI prediction may be 
misleading. 
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Table 1 : Statistics of the global fits for turbine design problem 




RS 2 

(Standard DSW) 

RS 3 

(Error-based DSW) 

RSquare Adj 

0.984 

0.998 

0.998 

rms-error Predictor 

151.947 

19.053 

22.103 

Mean 

-432.831 

513.476 

502.470 

Observations in Apay >0 / Observations, N 

4/9 

9/9 

9/9 

Testing rms-error 



133.616 i 


432 

432 

432 

Testing rms-error in Apay >0 

63.546 

22.645 

15.106 

#of the testing data in Apay >0, K 

193 

188 

188 

Mean of the testing data in Apay >0 

505.870 

563.886 

554.938 

Max. Error in Apay >0 

160.897 

45.356 

96.660 

Max. Error 

160.897 

885.188 

692.394 
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Table 2: RSI error at RS 2 data points for turbine design problem. Bold face rows are design points selected by 


Standard DSW. Average error of the boldface rows is 66.4 


D 

RPM 

Apay 



Error 

-0.5 

1 

212.66 

62.51 



0 

0 

15.75 

4.58 

0.497 

11.17 

0 

1 

806.18 

677.31 

0.832 

128.87 

0.3 

1 

983.04 

928.90 

0.792 




18.71 

16.22 

0.654 


i 

-0.6 

130.77 




i 

0 

630.5 



64.53 

i 

0.2 

749.29 

714.34 


34.95 

i 

1 

1074.38 

1173.87 


99.49 


Table 3: RSI error at RS 3 data points for turbine design problem. Bold face row are design points selected by 


D 

RPM 


)^ 51 


Error 

rn 

© 

1 



33.52 


74.89 

0 

0 


4.58 


11.17 

0 

1 


677.31 


128.87 

0.3 

-0.3 


14.34 

0.568 

9.74 

0.5 



703.76 

0.645 

25.13 

0.7 

1 


1127.53 

0.613 


0.9 

-0.6 

84.61 

24.46 


60.15 1 

1 

0 


565.97 

0.832 

ESI 

1 

1 


1173.87 

0.762 

99.49 


Table 4: Statistics of the RS fits forquartic polynomial 



RS 1 

(FCCD) 

RS 2 
(Standard 
DSW) 

RS 3 

(Error-based 

DSW) 

RS 4 

(RS 1 + RS 3 
design points) 

RSquare Adj 

1.000 

0.990 

0.847 

mmmMm 

rms-error Predictor 

0.000 

0.015 

0.056 

| 

Mean 

0.5800 

75164 

0.829 

0.653 

Observations in y >0.7 / Observations, N 

5/9 

7/9 

7/9 

7/9 

Testing rms-error 

0.08 

0.260 

0.259 

0.075 

#of the testing data, (M-N) 

432 

432 

432 

428 

Testing rms-error in y >0.7 

0.051 

0.073 

0.063 

0.039 

#of the testing data in y >0.7, K 

200 

198 

198 

198 

Mean of the testing data in y >0.7 

0.789 

0.794 

0.765 

0.776 

Max Error in y>0.7 

0.161 

0.149 

0.159 


Max. Error 

0.192 

0.778 

0.750 

0.217 
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Table 5: RS 1 error at RS 2 design points for quartic polynomial. Bold face rows are design points selected by 


Standard DSW. Average error of the boldface rows is 0.0916 



Table 6: RS 1 error at RS 3 design points for quartic polynomial. . Bold face rows are design points selected by 
error-based DSW. Average error of the boldface rows is 0.0868 



Table 7: Statistics of the RS fits for quartic polynomial in upper-right quadrant where testing rms-error in y >0.7 is 

0.061 by RSI 




RSquare Adi 


rms-error Predictor 


Mean 


Observations in y >0.7 / Observations, N 


Testing rms-error 


Testing rms-error in y >0.7 

0.026 

#of the testing data in y >0.7, K 

74 

Mean of the testing data in y >0.7 

0.790 

Max Error in y>0.7 

0.060 



For the comparison, test points are the designs of the relevant quadrant. 
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Table 8: RS 1 error at RS 5 design points for quartic polynomial in upper-right quadrant. Bold face rows selected 
by Standard DSW. Average error of the boldface rows is 0.0943 


0.7417 


0.7029 


.5975 

.7434 

.5621 


74 

75 


0.9380 


1.0000 


0.6615 

0.8315 


0.7148 

0.7617 


Error 


0.0000 


0.0798 


0.0003 


0.0001 


0.0003 


Table 9: RS 1 error at RS 6 design points for quartic polynomial in upper-right quadrant. Bold face rows selected 
bv error-based DSW. Average error of the boldface rows is 0.0658 


0.7417 


0.7029 


RSI 

E 

0.742 

0j 

0.703 

0. 


1.0000 


■ 

.6128 

— 

0.1781 

E 

.8315 

mm 

0.0003 

0.6128 

0.977 


0.7617 

1.000 



(a)FCCD (b) y A c con tours 

Figure 1: Design space for two-variable examples: FCCD design points and J A G con tours (Mean yj X Q =0.0646) 
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(a) Standard DSW approach (b) Error-based DSW approach 

Figure 2: Standard and error-based design space windowing approach for searching regions of high response designs 
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(b) Absolute error of RSI 


Figure 3: Apay contours and RSI errors for turbine design problem (x^D and x 2 : RPM) 
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(a) RS 2 (Standard DSW) 


(b) RS 3 (Error-based DSW) 


Figure 4: Design space for RS models for turbine design problem (x } :D and x 2 : RPM) 




(a) RS 2 (Standard DSW) (b) RS 3 (Error-based DSW) 

Figure 5: Error contours after DSW in turbine design problem (xj : D and x 2 : RPM) 




(a) RS 2 (Standard DSW) (b) RS 3 (Error-based DSW) 

Figure 7: Design space of global RS models for quartic polynomial 
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(a) Data points (b) Error contours 

Figure 9: RS4 for quartic polynomial (RS 1+RS 3 design points) 


(a) RS 5 (Standard DSW) (b) RS 6 (Error-based DSW) 

Figure 10: Design space of local RS models for quartic polynomial in upper-right quadrant 
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(a) RS 5 (Standard DSW) (b) RS 6 (Error-based DSW) 

Figure 1 1 : Error contours for quartic polynomial in upper-right quadrant 
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APPENDIX 


Response Surface Methodology 


Response surface approximations fit numerical or physical experimental data with an analytical model that is 
usually a low-order polynomial. The response at design point x* is denoted as y [ - y(x) ] and given by Eq. (4) 
y{x) = 7}(x) + e (4) 

where rj(x) is the true mean response at design point x and € represents other random sources of variation such as 
measurement error and noise not accounted for in rj . In other words if the experimenter has the luxury of 
performing the experiment many times at x, the average of the observations will tend to 7 } despite the random errors 
€ . Random errors £ are assumed to be uncorrelated and normally distributed random variable with zero mean and 
standard deviation <J , which is the same at all points. In RS technique the true mean response is assumed to be 
given in terms of coefficients /?, s and shape functions /;(x) (/=!, fib ) as in Eq. (5). 

**)=£fl/i00» f T P 

i=l 

where f T = {/,(x ) / 2 (x) f„ h (x )} and P = /? 2 ■■■ P „ h }' and T denotes transpose for the matrix 

representation. The n b shape functions f i are usually monomials and /?, are unknown coefficients, representing 
the best approximation toy when noise is absent. With noise the fitted approximation is given as 

X*) = & i f i (x) = f T b < 6 > 

1=1 

where bjS (vector b) are estimates of the P i (vector p ) obtained from a least squares fit. The difference (residual) 
between the datay for a point x and the estimate defined in Eq. (6) is given as 

e = y-y(.*) 

The residual for the N data points can now be written in matrix form, 


e = y-X 1 b 1 


( 8 ) 


where Xj is the matrix whose terms in the row associated with the point x are formed by monomials /, ( X ) *. For 
instance, for a quadratic model in two-variables with N data points X , is given as 


(9) 



1 

*n 

*21 

*n 

*11*21 

*21 


1 

*12 

*22 

*1 2 2 

*12*22 

*22 

x, = 

1 



2 


2 


x y 

X 2J 

*1, 

*1 y*2 j 

X 2J 


1 

* IV 

*2 N 

*IV 

X \N X 2N 

*2 N 

The coefficient vector b, 

in Eq. (8) 

is solved for minimum 

resi 


expressed as 

b, =(X7x,) , x7y 

When the model is exact (no bias error), unbiased estimates of the noise <7 in the data is given as 
_ ) y T y-b7x?y 

V N ~ n b 

With finite number of data, errors in the data cause errors in the coefficients, and that in turn causes a prediction 
error of the RS approximation that depends on location. Because Eq. (5) is only an approximation to the true mean 


( 10 ) 


( 11 ) 


* Bold face is used for vector and matrix representations, e g. x = [a , x 2 ]' for two-variable case where superscript 
T stands for transpose operation 

* The subscript “1” used for p u b , and X, refers to the set of monomials included in the model. Later, a subscript 
“2” will be used for the set of monomials needed to complete the model in order to obtain the true response. 
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response function, 5 will contain not only noise error, but also error due to the approximation. Besides .v , the 
quality of the approximation is often measured by the adjusted coefficient of multiple determination R 2 . 

2 


R: = 1 — 


( 12 ) 


T.(yj-y) 2 /W-i) 

7=1 

where y is the average value of the response data. 


D-Optimal Design 

A D-optimal design minimizes the generalized variance of the estimates, which is equivalent to maximizing the 
determinant of the moment matrix, M [15]. 

X T X 

v w / 

N n » 

The D-optimal design approach makes use of the knowledge of the properties of polynomial model in selecting the 
design points. This criterion tends to emphasize the terms of the polynomial model with the highest sensitivity [17]. 



Mean Squared Error Criterion [9] 

As a measure of the error in the approximation, the mean squared error of prediction MSEP defined as in Eq. 
(14) is used 

MS££(x) = £[P(x)-r 7 (x)] 2 (14) 

where 7/(x) and j)(x) are the true mean response and the prediction by the fitted model, respectively at x. MSEP is 
by definition an expected value that would be reached if the number of data used in approximation were unlimited. 
Equation (14) may be rewritten as 

MSEP = E\y(x) - Ey(x)] 2 +[£>-(x)-^(x)] 2 (15) 

The first term in Eq. (15) represents the variance error due to random noise and the second term represents bias 
error due to inadequate modeling. This error expression is usually averaged via integration over the design space 
and the integral is minimized by choosing the experimental designs (DOE) that control the effect of one or both 
types of error-noise and bias-[18 and 5]. It is reported in [18], “the averaging MSEP over the design region in fact 
mask a poor performance by the RS approximation at certain locations of the design region ” Papila and Haftka [9] 
instead attempted point-wise characterization of the error and to determine the design regions where RS prediction 
may suffer due to either or both types of error. Therefore Eq. (1 5) is used to investigate the variation of MSEP from 
point-to-point. 

The expectation of the predicted response at a given design point x can be expressed as 
Ey(x) = f?E(b t ) (16) 

where fj * is the vector of shape functions f i [see Eq. (5)] calculated at point x. 

The mean and the variation for the coefficient estimates are given as 

£(b 1 ) = (xj'x 1 }‘x7£(y) (17) 

Kar(b,) = a 1 (xj Xj f (18) 

where a is the standard deviation of the noise error. The first part of the mean squared error in Eq. (15) is equal to 
the prediction variance at x. 

£[y(x)-£j)(x)p =Var[y(x)] 

= f t T Var(b t )f| (19) 

= cr 2 f 1 T (x[X 1 )' 1 f 1 

The second part of Eq. (15) is the squared error due to the inadequacy of the model used 


* Subscripts 1 and 2 assign matrices and vectors for the fitting model and the missing terms, respectively. 
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[£p(x) “ 7(*)] 2 = (#(asLv’(x)]) 2 ( 20 ) 

Therefore mean squared error of prediction in Eq. (15) can be rewritten as 
MSEP = Kar[y(x)] + (fl/as[y(x)]) 2 ( 21 ) 

The form of the prediction variance, Var\y{x)\, in Eq (19) does not change even in the presence of inadequate 
modeling, since variance is only caused by random noise. However, when bias error is present, the residual error 
mean square 5 2 is not an unbiased estimate of <7 2 . Using Eqs. (10) and (1 1) we get 

s 2 _ y T ( 1 N ~ p )y (22) 

N-n b 


where P = X,(x?'X 1 ) _1 X?' ( 22 ) 

The accuracy of this estimate depends on factors such as the available data or number of points and also the 
adequacy of the fitting model that determines whether the estimate is unbiased or biased. 

The true mean response at x can be written as 

^(x) = f 1 T p l+ f 2 T P 2 ( 24 ) 

where f 2 are terms missing from the assumed model. Since one usually does not know the true response, it is often 
assumed to be a higher order polynomial when monomials are used as shape functions. For instance, for a quadratic 
fitting model in two-variables when true function is cubic 
f, T =[l Jtj x 2 x 2 x t x 2 X 2 ] 

-\Pu Pn Pn P\ 4 Pis P\6 ] (25) 

^7 = [ r l 3 x ^ x 2 X \ X 2 * 2 ] 

P 2 = \Pl\ Pll Pi 2 Pi 2 ] 

The mean of the true response is given as 


£(y) = X l p l +X 2 p 2 (26) 

where X 2 is similar to X 1? but due to the terms (shape functions) present in the true response that are not included 
in the fitting model. After substitution of Eq. (26), Eq. (17) becomes 

£(bj) = Pi + Ap 2 (27) 

where A = (X^X 1 )" 1 X^X 2 is called the alias matrix. Substitution of Eq. (27) into Eq. (16) yields 

£j)(x) = f 1 T (p 1 +Ap 2 ) (28) 

(^[Kx)]) 2 = |?i T (Pi+AP 2 )-f, T Pi -f 2 T P 2 J (29) 

= Pl(A T f,-f 2 ][f 1 T A-f 2 T |P 2 

The MSEP at a given point can now be estimated as by using s 2 instead of C 2 

MSEP(x) = s 2 f^(Xjx 1 y l U +PjMp 2 (30) 

where M^A 7 ^ -f 2 ]If, T A -f 2 T | ( 31 ) 


Note that MSEP is an expectation by definition [Eq. (14)], and Eq. (30) is its estimate by the data available and 
associated s 2 . Using £(y) from Eq. (26), the expected value of biased error mean square s 2 given by Seber [19] 


E{S 2 )=a 2 + m<hL^Mi 

N- Pi 


(32) 


since o N -p) is an idempotent matrix, that is a matrix whose square is equal to itself. It is also positive semi- 
definite matrix, so E(s 1 )'2.o 2 provided that (I N -P)X 2 p 2 *0. Equation (32) means that when bias errors are 
present, s 2 is a biased estimate of O 2 . 

If this expected value, E(s 2 ), is substituted in prediction variance contribution in Eq. (30) expected value for the 
estimate MSEP can be expressed as 
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E[MSEP(x)] = f?(X]X t y'fi +pjMp 2 

N-p x J (33) 

so-^cx^x.r'f.+pjGp, 

where K = X ] (X 2 - X, A) (34) 

G= »Ak + M (35) 

N- Pi 

The aim of the method is to identify points where the bias error is large. It is assumed that the terms missing 
from the fitting model are known, but there is not enough data to calculate the corresponding coefficients P 2 ■ If one 
can estimate the size of JJ 2 , it is possible to formulate a constrained maximization problem for the largest magnitude 
of the mean squared error predictor that may be experienced at any given design point for the worst possible p 2 °f 
that magnitude. 

max£[A/S£P(x)] 

(36) 

such that ||p 2 1 —C 

The Lagrangian for this optimization problem can be written as 
Z.(P 2 ,A) = £[M§£P(x)] + ,l(pJp 2 -c) (37) 

Differentiating the Lagrangian with respect to P 2 

V|<r 2 f, T (X?'x 1 )- , f 1 l + V(pjGp 2 ) + AV(pJp 2 — c) = 0 (38) 

T a a 3 f 

where V = • * • . Equation (38) yields the following eigenvalue problem at a design point 


Gp 2 + Ap 2 = 0 or 

CP 2 -A G p 2 =0 (39) 

for which the maximum eigenvalue ^ max characterizes the maximum possible mean squared and bias error 

associated with the assumed true model that includes shape functions missing in the fitting model. The 
corresponding eigenvector defines the coefficients of the missing shape functions that results in the largest bias error 
when fitted only with the assumed model. The eigenvectors and the function experiencing the worst possible bias 
error may be different point-to-point although the magnitude of the missing coefficient vector is constrained. So the 
eigenvalue calculated does not reflect the true function corresponding to the data (as the data is insufficient to 
calculate P 2 ). It reflects instead, the assumed form of true function with the shape function coefficients p 2 (among 

all the possible combinations such that pjp 2 = c ) causing the largest error. 

As can be seen from the derivation the eigenvalue error measure strongly depends on the DOE used, but not on 
the response data. This independence from the response data may be misleading especially for the cases where 
fitting model predicts the true function very well at the data points. The following 2D cubic polynomial example 
demonstrates such a case. 

y = \+xt+ X l (40) 

Face-centered central composite design (FCCD) in 2D was used for constructing a quadratic RS. Figure 1 
presents the FCCD points and relevant yjZ G field. The maximum values are at the data points on the boundary, but 

Table A1 indicates a perfect with the nine data points. The testing rms-error (using 21 by 21 grid points except the 
nine data points) and the maximum error, however, shows that the fit is actually quite poor. Figure A1 presents the 
true function and the error field. Comparing Figure Alb and Figure lb indicates that eigenvalue error measure did 
not predict the high-error regions in this particular example, unlike for the examples studied and reported in [9]. It is 
expected that the benefit of the error measure will be increased further if the response data can also be used in the 
derivation. 
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Table A1 : Statistics of the quadratic RS for cubic polynomial given in Eq. (40). RS is given is y - 1 + x l + x 2 


RSquare 

1.000 

RSquare Adj 

1.000 

rms-error Predictor 

0.000 

Mean 

1.000 

Observations, N 

9 

Testing rms-error (21 x 21 grid) 

0.385 

Max. Error 

0.768 




(b) Absolute error of RS 1 
for cubic problem in Eq. (40) 
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